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Abstract 

What input signals will lead to synchrony vs. desynchrony in a group of biological oscil- 
lators? This question connects with both classical dynamical systems analyses of entrainment 
and phase locking and with emerging studies of stimulation patterns for controlling neural 
network activity. Here, we focus on the response of a population of uncoupled, elliptically 
bursting neurons to a common pulsatile input. We extend a phase reduction from the lit- 
erature to capture inputs of varied strength, leading to a circle map with discontinuities of 
various orders. In a combined analytical and numerical approach, we apply our results to both 
a normal form model for elliptic bursting and to a biophysically-based neuron model from the 
basal ganglia. We find that, depending on the period and amplitude of inputs, the response 
can either appear chaotic (with provably positive Lyaponov exponent for the associated circle 
maps), or periodic with a broad range of phase-locked periods. Throughout, we discuss the 
critical underlying mechanisms, including slow-passage effects through Hopf bifurcation, the 
role and origin of discontinuities, and the impact of noise. 

Keywords: elliptic bursting, circle maps, perturbed oscillators, synchrony, mathematical neuro- 
science 

AMS subject classification: 92B25, 92B20, 34C28 

1 Introduction 

Many types of physical and biological systems exhibit intrinsic bursting - rapid discharges of con- 
secutive, fast dynamical events separated by periods of quiescence. In particular, bursting neurons 
serve myriad functions in the nervous system; prominent among these is their role in central pattern 
generators that create rhythmic neural activity [131 UHl EQ US] . Bursting dynamics also feature in 
pathological oscillations associated with disease conditions, as for basal ganglia networks in Parkin- 
sons disease [56j ESI til EH] , where elevated synchrony and rhythmicity among neurons is linked 
to motor symptoms. 

Here, we focus on synchrony and desynchrony among bursting neurons in the simplest possible 
setting: a population of uncoupled bursting neurons receiving a common input signal. We study 
elliptic bursters [38j [57] - that is, non-linear oscillators with fast and slow variables, and for which 
burst onset is caused by passage through a (subcritical) Hopf bifurcation in the fast subsystem and 
burst offset follows from a saddle-node bifurcation of limit cycles (see Sect. [2] below). The driving 
signals are periodic pulsatile inputs. 
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We find that there is rich variety in the response to these inputs, depending on their strength 
and frequency. An illustrative example is presented in Fig. [TJ where we plot simulated voltage 
traces of two bursting cells, both receiving a common pulsatile input I(t). In the left panel, the 
cells' bursting phases are initially well separated, but a pair of "strong" input pulses synchronizes 
them. In the right panel, the cells' phases are initially nearly synchronized, but a relatively "weak" 
input drives them apart. As we will show, the outcome depends on pulse strength and on inter-pulse 
timing in interesting ways that arise directly from the dynamics of elliptic bursting. 
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Figure 1: Simulated voltage traces of two uncoupled cells receiving common inputs (from the 
conductance-based neuron model of Eqn. (38); see Sect. [5]). (a) A strong input synchronizes cells 



that are initially out of phase, (b) A weak input desynchronizes cells with close by initial phases. 



Our goal is to understand the mechanisms responsible for this and related phenomena. Specif- 
ically, we ask whether or not cells will entrain or become separated under the driving effect of a 
common periodic input signal. Our main goal is to develop and explain a general answer to this 
question. We note two broad scientific motivations for doing so, although exploring the implica- 
tions of our results for each largely remains for future work. First, entrainment of a population 
of uncoupled cells to an input signal determines the "reliability" of a neuron's response - that is, 
the repeatability of a response to a fixed input signal on multiple trials in which the neuron is in a 
different initial condition. Reliability is fundamental in understanding how neural dynamics encode, 
e.g., sensory signals [HJ Ell ESJ [55j E0, [4] . Second, common input signals to populations of bursting 
cells (of elliptic and other types) occur naturally in layered neural networks, as in the basal ganglia; 
in this brain area, common, pulsatile electrical stimuli are also artificially applied as a therapy for 
Parkinson's disease j56j [66l [7j. We give a very brief application to this general setting in Sect. [6| 

Throughout this paper, we use a normal form model for elliptic bursting developed by Izhikevich 
[38] . It is the simplest system that captures the dynamical features of elliptic bursting and we show 
that it accurately describes a more complex bursting neuron model derived from basal ganglia 
physiology. Thus, we will often refer to the normal form model as describing a "cell." 

Building from, e.g. [7J, we show that the dynamics of the normal form model under periodic, 
pulsatile inputs admits an accurate reduction to a discontinuous circle map that can be analytically 
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defined. This map forms the basis for our theoretical results, and correctly predicts synchrony 
and desynchrony for the physiological model. Nevertheless, because our results are linked to the 
normal form model, we note that our approach and results have potential applications in the study 
of general slow/fast oscillators undergoing a delayed bifurcation, well beyond neuroscience. 

The paper proceeds as follows. Section [2] deals with the analysis of elliptic bursting dynamics as 
well as phase reduction to a discrete dynamical system on the circle. Section [3] presents an analysis of 
the reduced dynamics and explores synchronizing and desynchronizing effects of common pulsatile 
inputs. Next, in Section [4j we study the effect of noise, which has a non-trivial and interesting 
impact on the circle map and resulting patterns of synchrony. In Section [5| we carry out a series 
of numerical experiments that verify our reduced models. Finally, in Section |HJ we show how the 
circle map framework can be used to analyze the effect of multiple sequences of pulsatile inputs. As 
a proof of concept, we present a brief example in which an input signal is designed to compete with 
the effect of a global entraining drive synchronizing a population, as in the basal ganglia setting 
described above. Our principal finding - unifying analyses of several models and settings is that a 
population of (elliptic) bursting cells will desynchronize in the presence of weak to moderately strong 
common inputs, if these inputs have a frequency slightly slower than the natural burst frequency. 

2 Geometry of forced elliptic bursters 

In order to better understand the effect of a common stimulus on a population of bursters, we 
first describe the dynamics of single cells. Dynamical models that capture bursting usually include 
multiple timescales and are often called slow/fast systems. Indeed, most intrinsically bursting 
solutions arise from the evolution of one or more slow variables that periodically steer fast variables 
into distinct dynamical regimes - here, spiking and resting. 

2.1 Timescale dissection and basic model 

Slow/fast systems can be written in the form 

z = f(z,y,e) 
y = eg{z ) y,e) 

where z is a vector of fast variables, y is a vector of slow variables and e is the slow/fast timescale 
ratio. Such systems arise in many areas of mathematical modeling and can describe general multiple 
timescale phenomena. In the singular limit, where e —> 0, one obtains an equation where the slow 
variable(s) y can be considered as parameter (s) for the fast subsystem (z). This approach allows 
one to investigate the dynamics of the fast subsystem and subsequently to construct solutions of 
the full system by carefully reintroducing slow dynamics - i.e., by studying the perturbed system 
(e / 0). An elegant mathematical toolset known as geometric singular perturbation theory has 
been developed to study such phenomena [T71 [43] . 

The results brought forward in this paper relate to the effect of pulsatile perturbations on 
slow/fast bursters in which a delayed Hopf bifurcation is central to the onset of rapidly varying 
dynamics. One of the first systems of this type to be studied was a bursting phenomenon in the 
Belousov-Zhabotinskii chemical reaction [59]. The mechanisms central to our study can be found 
in many chemical, physical and biological systems. 
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However, our main purpose is to better understand the effect of perturbations on conductance 
based models of single neurons (i.e., models of Hodgkin-Huxley type) that possess such a separation 
of timescales. In many cases, calcium concentration acts as a slow variable while voltage and 
associated ionic currents evolve on the fast timescale. Rinzel and Lee first studied fast/slow solutions 
to models of parabolic bursters in 1987 using singular perturbation methods [58J. Since then, much 
effort has been invested in understanding bursting solutions arising in conductance-based neural 
models and their reduced forms [EH [57] . In 2000, Izhikevich produced a classification of bursting 
mechanisms [36], including all of the possible codimension one bifurcations of the fast subsystem 
that could be responsible for the onset and termination of spiking dynamics. 
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Figure 2: Elliptic bursting trajectory from model ([2]). (a) Trace of Re(z). (b) Trace of y. (c) 
Projected into the {Re{z) ) y) plane, thicker black lines show y-parametrized fixed points and limit 
cycles of the fast subsystem (i.e., from the singular limit). S is a family of fixed points (stable 
for the solid line and unstable for dashed), U is a family of unstable orbits and P is a family of 
stable orbits. Thin black line shows bursting trajectory, with red arrows indicating time evolution 
direction. 



We concentrate on cells characterized as (subcritical) elliptic bursters or type III bursters, in 
which a subcritical Hopf bifurcation of the fast subsystem drives the onset of bursting and a saddle 
node bifurcation of limit cycles is responsible for the return to silent state (see Fig. [2]). This 
type of intrinsic bursting is well understood in the absence of forcing (input). In particular, Su, 
Rubin and Terman established existence and stability properties of elliptic bursting solutions in [63] . 
Izhikevich [38] presented a portrait of elliptic bursting dynamics by describing the fast subsystem via 
the normal form of the (codimension two) Bautin bifurcation, and derived a closely-related model 
of weakly coupled networks of these bursters [39]. Throughout this paper, we keep our analysis 
as general as possible and often illustrate our results with a variation of Izhikevich's normal form 
model (see Eqn.Q) to simplify mathematical manipulations. However, we stress that our analysis 
can be carried out for any elliptically bursting model. 
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The normal from model is 

z = (y + iw)z + 2z\z\ 2 - z\z\ 4 + I(t) 
y = e(a- \z\ 2 - by) 

where z G C represents the fast variable and i/GR the slow variable. This is as found in [39J, but 
with the term by added to the slow variable dynamics in order to explore the effect or nonlinear 
ramping in the delayed Hopf bifurcation, a feature that is found in many models for which bursting 
dynamics are caused by slowly varying calcium concentration. We set w = 1 and e = 0.01 for 
the remainder of the paper and will consider distinct cases in which we vary a and b. Intrinsically 
bursting solutions generally arise for parameter choices where y is positive when z is at rest and 
negative when z is spiking (oscillating). 

Notice the forcing via the signal I(t) in the equation for the fast variable. We wish to model an 
input signal that causes an instantaneous voltage response, so we set I(t) = ^2 n A n S(t — t n ) where 
5 is a delta function and A n G M + . We will refer to these perturbations as t n - kicks of amplitude 
A n \ a kick simply translates a solution at time t n by an amount A n in the real "direction" of z. In 
this paper, we focus on periodic kicks of fixed amplitude A = A n and equal spacing r = t n — t n _\ 
between kick times. 



2.2 The elliptic bursting cycle 

Figure [2] shows an elliptic bursting trajectory from numerically integrating Eqn. ^ with a = 0.8 
and 6 = 0. As mentioned above, the standard approach is to think of the slow variable y as a 
parameter that determines the fast dynamics (z). The fast subsystem undergoes a subcritical Hopf 
bifurcation at yn — and a saddle-node of limit cycle bifurcations at y$N = — 1- When e — 0, and 
given a fixed y G (ysN, Vh): the fast subsystem is bistable and has a sink at z = inside an unstable 
periodic orbit, itself contained within a stable periodic orbit. These y-parametrized limit sets form 
normally hyperbolic invariant manifolds that persist under small perturbations (1 >> s > 0) [63] . 
We denote by S the family of equilibria (z — 0). These are stable when y lies to the left of yn 
and unstable on the right, together forming the silent branch. P is the family of stable periodic 
orbits and represents the spiking branch. For Eqn. ([2]), the radii of these orbits about S are given 
by rp(y) = y/l + y/y + 1 for y > ysN- Finally, U refers to the family of unstable orbits with radii 
ru(y) — a/1 — Vv + 1 f° r y ^ (vsn, Vh), acting as a separatrix between the stable side of S and P. 
A bursting solution occurs when the slow dynamics of y steer the fast subsystem rightward along 
the silent branch S ) until the Hopf point y# is reached. The solution then keeps moving rightward, 
sticking close to S for some transient period even though the equilibria forming S are now unstable, 
but is eventually attracted to P when y reaches yj, where spiking begins. The y dynamics then pull 
the oscillating fast subsystem leftward along P, until the latter vanishes at ysN, where the solution 
is attracted back to S and another cycle begins. 

Of particular interest is the slow passage effect through the Hopf point in which the solution 
does not immediately jump up into spiking when S looses stability (yj ^ yn)- This delayed 
bifurcation phenomenon as been previously studied [3j El [63] and its implications in the context of 



pulsatile perturbations will be established in Sect. 2.4, Importantly, several authors have shown that 



noise can sharply diminish this slow passage effect [48j [471 El E3] . While we first treat the noiseless 
case, we study the effect of stochastic terms on the phase reduction and response dynamics in 
Sect. [1 
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2.3 Phase reduction 



There is a substantial body of literature concerning phase reduction of oscillators and their behavior 
under noise, forcing, or coupling [24J, [4Ql HOl UHl [30l [53], [65J, [23] . The general idea is simply to associate 
the endpoints of a periodic solution's cycle and to parametrize the movement along this solution 
with a phase 9 e S 1 . Although the point 9 = is arbitrary, it is often chosen to correspond to a 
distinguishable event within the periodic orbit, such as the apex of an action potential in a model of 
a spiking cell. This reduction becomes useful when the limit cycle has some stability properties and 
one can track the phase response of the solution following a perturbation: specifically, by computing 
the phase difference on S 1 between the unperturbed solution and the perturbed one, as t — » oc and 
the latter contracts back to the limit cycle. When well defined, this one-dimensional description has 
the advantage of being analytically tractable while preserving the behavior of an oscillator subject 
to perturbations. 

For systems with asymptotically stable limit cycles, phase reduction can be carried out rigorously 
as long as the effect of forcing translates the solution to a point in the cycle's basin of attraction. 
This basin is foliated by the strong stable manifolds of each point on the limit cycle. By knowing 
on which manifold - or isochron - the solution lands following a perturbation, we know exactly to 
which phase it will be attracted in the limit as t —> oc [241 EH EH E21 dU |37j . 

We next describe a phase reduction for elliptic bursters. As we will discuss further, the bursting 
trajectories are not stable limit cycles, and so we cannot directly compute isochrons. Nevertheless, 
the difference between the timescale of the burst period and the timescale of attraction normal to 
the (singular limit) solution enables us to proceed. We closely follow Best et al. [7J, who derive 
circle dynamics for an elliptic burster receiving periodic inputs from a model excitatory neuron. In 
their work, each excitatory kick always transitions a solution from the resting to the spiking state 
- i.e., to the branch of periodic orbits P in Fig. [2| if it is not already following that branch. Best 
et al. also use an approximation of linearity for the slow trajectories (i.e., y is piecewise constant). 
Here, we relax both of these assumptions, in particular while studying the response to weaker kicks 
that do not necessarily generate a burst. As we will see, it is these weaker kicks that will lead to 
desynchrony; that is, dynamics that appear chaotic or are phase locked at high period. 

In [63], the authors use Fenichel theory to show that there exist 0(e) neighborhoods N$ and 
Np around S and P such that - if a solution enters either the left side of Ns or the right side of 
Np and the slow dynamics behave as mentioned above - then the solution will transition between 
the two neighborhoods in a periodic, bursting fashion. Furthermore, they use averaging techniques 
to show that the period of such a cycle can be approximated up to 0(e) by the sum of the passage 
times Ts and T P along the respective manifolds. 

Although it is not clear whether or not there exists a single periodic solution of the full system, 
[63] shows that such dynamics are at least metastable. That is, Ns and Np are locally attracting 
and any bursting solutions must live in these sets. Furthermore, solutions starting outside the two 
attracting sets are quickly attracted back to them. 

Numerically obtained solutions of Eqn. ^ do not trace back exactly the same path from 
cycle to cycle but have periods that vary minimally, as expected [63j. Specifically, we numerically 
integrated a solution in order to obtain 150 burst cycles and computed the coefficient of variation 
( Cv = standard^deviation ^ of the cyde durations> W e found CV = 0(W~ 3 ) for parameter choices 

yielding bursting solution in Eqn. ^ (numerical methods as in Sect. [5]). Whether there is a 
periodic solution with a much longer period than the bursting cycle, or whether solutions are 
instead quasiperiodic or aperiodic remains an open question. 
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The small CV (indicating a robust cycle period), along with the metastability described above, 
motivate an approximate reduction to dynamics on the circle, as in [7J. We will revisit the notion 
of uncertainty in cycle periods in Sect. [4j However, for what follows, we use the singular limit 
assumption that a bursting trajectory evolves along S and P with well defined passage times T# 
and Tp, and will use this trajectory to compute phase reduced dynamics. 




Figure 3: Schematic presentation of phase reduction. We associate endpoints of a bursting cycle 
(where y = ysN), and map the trajectory onto the unit circle. 

As illustrated in Fig. [3j we represent bursting solutions by a phase variable 6{t) G S 1 = R/Z. 
We let 9 = ~ 1 correspond to the "jump down point" on the bursting trajectory, where solutions 
transfer from spiking to resting. We choose this reference point because the spiking to resting 
transition is fast and is associated with a constant value of the slow variable, y = ysN (unlike, as we 
will see, the transition to spiking following a pulsatile input). Essentially, the phase 9 of a bursting 
cycle is given by time rescaled by the period (9 = Tg + Tp ) where at t = 0, y = ysN- Note that 
the first portion of the unit circle following 9 = represents the silent branch S and the remaining 
portion represents the spiking branch P. 

2.4 The kick map 

We now study phase dynamics of bursters receiving pulsatile inputs (kicks) . Specifically, we derive 
a phase translation mapping Fa{9) ) such that if 9 is the phase of a cell when a kick of strength 
A arrives, Fa(9) is the phase of the kicked solution - relative to the unperturbed solution - after 
it relaxes back to the burst cycle. We will refer to this as the kick map. In [7J, a similar map 
for elliptic bursters is derived, and this idea inspired the present work. However, there are two 
differences with the map we derive here. First, the phase of trajectories in [7] is defined relative to 
the period of pulsatile inputs; in our case, phase is defined relative to the (unperturbed) period of a 
burst trajectory. This latter construction has been previously used in the context of integrate and 
fire cells with soft reset [5] , and we find that this makes it easier to visualize the role of changing kick 
amplitude and period on the structure of the map. Second, in [7], only "strong" kicks are considered; 
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as we will see, the map develops additional features, including discontinuity and expansion, in the 
case of weaker kicks. 

As discussed in Sect. |2.3[ we assume that unperturbed elliptic bursting solutions have fixed 
times spent in silent (T#) and spiking phases (T P ). When computing a map for a given system as 



done in Sect. 2.5, it is best to work with unsealed time and later (implicitly) normalize the time 
variable by the burst period (T# + T P ) so that our phase variable 9 remains between zero and one. 
In this section however, we derive the kick map for an arbitrary elliptic bursting model and assume 
that the period is already unitary (T$ + Tp — 1) for simplicity. The rest of the notation follows 
that of system ^ but the reader should keep in mind that z and y represent general fast and slow 
variables. 

Our computations are intimately linked to the evolution of the slow variable y along the branches 
S and P. Recall that y spans [ysN,yj] when z travels along S or P. To better track the variable 
2/, we label its dynamics along S by y(t) = hs(t) and along P by y(t) = hp{t). 

Thus, for a burst trajectory that starts at y = ysN when t = 0, we have 

where h$(0) = ysN = hp(l) and hs(Ts) = yj = hp(Ts). Here, hs and hp are functions with ranges 
[VsniVj] an d respective domains [0,T#] and [T#, 1]. We assume hs is an increasing function while 
hp is decreasing. We now define the phase 9 of a bursting solution (z(i), ?/(i)) by 

nf h- s l {y{t)) if silent (z(t) e S) , , 

\ hp\y{t)) if spiking (*(*) eP). [V 

For unperturbed solutions, 9 — t mod 1. In general, expressions for hs and hp can be hard to find. 
As in standard approaches, one can integrate the dynamics of y by restricting the fast variable to the 
invariant manifolds S and P, and using the averaged motion of z on those manifolds f57l l23l l2l [63]. 
These calculations yield explicit formulas for the normal form system ([2]), as we demonstrate in 
Sect. ESI 

We are now equipped to define the kick map. Our first assumption is that a kick that arrives 
during the spiking regime has no effect on 9{t). Due to the separation of timescales, trajectories 
are attracted back to the stable limit cycles that form P in a vanishingly short time (with respect 
to the timescale over which the phase evolves). On the other hand, a kick while the cell is silent 
(z near S) may have distinct outcomes. If y(t) G [yH,yj] (sticking to the unstable part of S), any 
kick will send the cell to the spiking state since the trajectory is highly sensitive to perturbations. 
However, if y(t) G [ysN, Vh]-> one of two things can happen. If the kick is strong enough to translate 



the solution past the separatrix J7, then the cell jumps to the spiking state (Fig. [4 
other hand, the kick is not strong enough, the solution will attract back to S (Fig. 



(b)). If on the 
3(a)). 



9 w (A) = h- s \y w (A)). (5) 



From now on, we refer to a kick as strong if A > rjj{y) for all y G [ysN, Vh] where, ru(y) is 
the distance between S and U at y, in the direction of the kick. In other words, a strong kick will 
immediately result in spiking independently of the cell's phase. In contrast, we define a weak kick 
as one with an amplitude A < rjj{y) for values of y in some subinterval of [ysN> so that the 
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Figure 4: Weakly kicked (A = 0.5) trajectories from §2§. Top to bottom: Re(z) trace, y trace and 
y — Re(z) solution curve (red) truncated at t — 200 for clarity. Star indicates a kick and dashed 
blue line the kick's amplitude A. (a) Kick received at t = 80 and does not clear separatrix. (b) 
Kick received at t = 110 and clears separatrix. 



kick does not always immediately cause a cell to spike. (Note that a weak kick will result in spiking 
for any value of y outside this interval.) 
For the strong kick case, the kick map is 



While in silent phase, 9 G [0,T#], y increases according to h$(t)] when spiking is induced by a 
kick, the value of y is left unchanged but its dynamics are "reversed" and it starts decreasing via 
hp(t). The cell will then spike until y reaches ysN, which takes less time since we started closer 
ysN- Thus the impact of the kick is to advance the phase. This explains the first line of Eqn. ([6]), 
and is sketched in the panel (a) of Fig. [5| As already discussed, the kick has no effect when the cell 
is spiking, as expressed in the second line of Eqn. ([6]). 

For weak kicks, the situation is more complex. Recall that the branches U and P meet at y$N 
(^u(vsn) = rp(ysN)), and U vanishes at yn {tu(vh) = 0); we assume that ru is a non-increasing, 
continuous function. As a result, for rjj{ysN) > A > 0, we can find y w (A) such that A = rjj{y w {A)) 
and A < rjj{y) for all y G [ysNiVw(A)]. If the cell is in silent phase, y w (A) is essentially a cutoff 
point before which a weak kick will not elicit immediate spiking, as illustrated in Fig. [i](a). A weak 
kick delivered at any other point through the burst cycle will result in instantaneous spiking, as in 
the strong kick case. We recast this condition in phase coordinates, obtaining the cutoff phase 

What happens to a cell's phase when a weak kick does not evoke spiking (9 G [0, 9 W (A)])1 The 
trajectory is attracted back toward S but jumps up into spiking before it reaches yj ) as if it retained 
a memory of this past weak kick (Fig. [4] (a)). To better understand this phenomenon, we first need 
expressions for slow passage times through Hopf points. We use results for delayed bifurcations 




(6) 
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derived in [3] in order to predict points of transition between silent and active states. 

Let X(y) be the extremal eigenvalue of the fast subsystem linearized about the equilibrium points 
z (y), for some chosen value ofy G [ysN, Vj\- (Recall that this collection of points forms the manifold 
S = {zo(y)\y G [ysN,yj]}-) Points on S to the left of yn are sinks with Re(X(y)) < 0, whereas 
to the right, they are sources and Re(X(y)) > 0. Assume ^^-\ y=yH ^ 0. As a solution is pulled 
to the right by the slow dynamics, y crosses Re{X{y)) changes sign and solutions switch from 
being attracted to being repelled by S. However, this repulsion is not immediately apparent: the 
difference in time scales of our slow/fast system induces a discrepancy in spatial scales. As y varies, 
z is attracted to S for y < y H and repelled for y > y H , both at an exponential rates (proportional 
to Re(X(y))) on the fast timescale. 

The length of the slow passage (also called delay) to the right of yu depends on the dynamics 
of y along S. Borrowing notation from [3], suppose we can write 

y(t) = yi + 9(et) (7) 

where we assume g is a non-decreasing function and g(0) = 0. In this context, hs(t) = ysN + g(^t). 
If the system's full solution starts at (2^, yi) such that yi G [ysN, Vh] an d z % is f ar enough from z (yi) 
but is still in its basin of attraction, then the jump up point yj can be implicitly computed via 

= J {-g-\y- yi ))Re{Ky))dy (8) 

where g~ x is the inverse of g. In the unperturbed case, the point yj can be derived using yi = ysN 
in (|8j); this is sometimes called the memory effect for elliptic bursters [58] • We refer the reader to 
the appendix of [3j for the derivation of this integral condition. 




Figure 5: Schematic representation of the kick map 6 —> Fa(0) on the unit circle, (a) Strong kick 
map ([6]). (b) Weak kick map of Eq. ([9]). The solid arrows represent instantaneous change of phase 
while the dotted arrow indicates the evolution of the phase along the circle in time. 

If a weak kick does not elicit instantaneous spiking, we find that the "memory" starts anew at 
the time of the kick; in other words, if the kick is administered when y G [ysN, yw(A)} : then we set 
yi = y in Eqn. ([8]), and denote the associated jump up point by yj(y). In phase coordinates, if the 
cell is kicked at 9 G [0, 6 W {A)] ) y is given by hs{0). The onset of spiking happens at yj = yj{hs{0)) ) 
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or equivalently at 9j = h s 1 (yj(hs(9))). We capture this via: 

r + hp\yj{h s {fi))) - hs'iyjihsie))) if e [0, 6 W (A)) 
F A (0) = { h p 1 oh s (0) if 0e [9 W (A),T S ] (9) 

(e xde[T s ,i\. 

The first conditional definition maps the slow variable to its jump up value yj = yj{hs{9)) ) and then 
applies the strong kick map hp 1 o hs(0), and finally translates back by the phase —h^ 1 (yj(hs(9))) 
to account for the time taken in "slow passage" from yi = hs{9) to yj = yj(hs(9)). This mapping is 
sketched in the right panel of Fig. [5} note that it is only valid if the cell does not receive additional 
kicks before it enters the spiking state. 

This construction implies that the shape of a kick map does not vary continuously with the 
strength of a kick. Moreover, the only criterion which dictates the qualitative shape of the map is 
the value of 9 W (A) (we sometimes drop the A and write 9 W ). Other than determining the threshold 
9 W) the perturbative role of A is dimensionless since a kick acts on the fast variable as opposed to 
the phase which is defined over the slow timescale. That is, if two kicks of distinct amplitudes yield 
the same strong (resp. weak) outcome, the discrepancies between the times needed to attract close 
to the unperturbed trajectories are negligible. We emphasize that 9 W decreases as A increases: for 
strong kicks with A > ru(ysN), 9 W is always zero and the map does not change shape as A increases 
further. We use Eqn. ^ as the general expression for our kick map. 

We stress that there is a fundamental difference between the maps induced by strong and weak 
kicks. A weak kick always induces expansion in the kick map, as long as it speeds passage (yj ^ yj) 
through the Hopf point. Indeed, notice that by construction, hp 1 (y) > h^ 1 (y) for any y £ (ysN, yj) 
and hp 1 (y) = h^ 1 (y) when y = ysN,yj> From expression ([9]), we see that Fa(0) > 9 on (0,6 W ) 
and Fa(0) = 0. It follows from the mean value theorem that > 1 on some region contained in 
[0, 9 W ). We also note that when hp and h$ have similar shapes, it generally implies expansion of Fa 
on the whole interval [0,9 W ). To better illustrate this and other features, we compute expressions 
of this map for system ([2]). 

2.5 Computing the kick map for the normal form model 

In this section, we derive an analytical approximation of the kick map for the elliptic bursting normal 
form model. The task at hand is simple: use Eqn. ^ to compute the ingredients of expression 
h s (6), 4" 1} (#), h P \0), y 3 {y) and W (A). 

We first turn to h$ and hp, which are essentially the y components of a solution to Eqn. 
([2]), in silent and spiking modes respectively. In contrast with the previous section, we carry out 
computations using unsealed time which implies a full burst period X 7^ 1. One can still think of t 
as 9 in what follows, but the expression of the final map has to be rescaled. 

We exploit the separation of timescales in our equation and make the assumption, as in the 
singular limit, that z evolves exactly on the manifolds S and P. Notice that the y dynamics depend 
linearly on \z\ 2 . By substituting \z\ 2 = for hs, and \z\ 2 = rp(y) for hp, we obtain two scalar 
O.D.E.s 

y = e{a — by) when z is on S (10) 
y = s(a — rp{y) 2 — by) when z is on P. (11) 
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Bursting occurs when the right hand sides of (10) and (11) are respectively positive and negative, 
steering the fast dynamics in the required directions along S and P. Here we concentrate on 
parameter values a > and 6 > which satisfy this condition. As we will shortly see, 6 > implies 
that the evolution of y along S follows a saturating exponential ramp while 6 = implies a linear 
ramp. Both scenarios are found in biological systems, and - as we now show - the resulting kick 
maps have common characteristics which unify their response to pulsatile perturbations. 



For Eqn. (10), we can easily solve and get 



h s (t) = eat + Cs 



6 = 
6 > 



where Cs is an integrating constant. Setting y = ysN(=— 1) at t 
for d — and Cs = Vsn — t for 6 > 0. In turn, we have 



(12) 
(13) 

0, it is easy to see that Cs = ysN 



hs\y) 



y - vsn 

ea 



6 = 



h s 1 {y) = H4 a - h}) + ^ H4 a - WD b > o. 



(14) 
(15) 



For Eqn. (11), recall that rp(y) = + ^fy + 1 when y > ysN- Solving this O.D.E. is not as 
straightforward; fortunately, we only need the inverse hp 1 (which is a proxy for t in this context) 
to compute our map. This can be obtained directly via integration: 



hp\y) = 



eb 



hp\y) = --[(a - 1) ln(-a + ^+1 + 1) + y/^+T] +C P 6 = 



2 tan 



26y^TT+l 



v /-4(a-l)6-46 2 -l 



y/-A{a - 1)6 - 46 2 - 1 



ln(-a + by + ^y + 1 + 1) 



C P 6 > 0. 



Using (14), (15) it is straightforward to compute the time (T$) it takes for h$(t) to reach yj, and 
then derive values for Cp such that hp 1 (yj) = T# as required by our definition of h p . 

To compute the jump up point yj, we derive an expression for y^Vi) which is also needed in the 
definition of our kick map. In order to use the integral condition ([8]), we must first write expressions 
for the y dynamics in the form of Eqn. Q. Here we substitute ysN by yi in the expression of Cs 
for (12) and (13) to allow for arbitrary initial conditions and get 



y(t) =yi + ast 6 = 
I/ (t)= y< + (l-e- tet )(|-yO 6>0 



which yields 



g(st) aet 



-bet 



)(' 



6 = 
b > 



(16) 
(17) 
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which in turn give 



g-\y) 



y 



a 



ln(l + — 

g-\y) = - r a/ b>o. 



(18) 
(19) 



It is easy to compute the extremal eigenvalue for the linearization about z = of the fast 
subsystem in Eqn. ^ : X = y ± iw and therefore Re(X(y)) = y. Turning now to the integral 
condition ([8]), when b = we use (18) and write 



Vi 



a 



-dy 



by which we can deduce that 



0. 



(20) 



In other words, when the slow variable follows a linear ramp along the silent branch S (see (16)), 
the jump up point yj is symmetric to the initial point yi about yn — 0. This can be seen from Fig. 
[2] where y S N -1 and yj 1. 



In the case where b > 0, y follows a saturating exponential ramp along S (see (17)) which implies 
that y decelerates as it moves rightwards and thus shortens the slow passage. Using (19) in the 
integral condition ([8]), we get 







y 



dy 



which gives 



thus implying the relation 



We then isolate yj to get 



1 a\n(by — a) 



+ y 



Vj 



y=yi 



b byi-a 



b>0 



(21) 



where W is the Lambert W function (product logarithm function). 

The last ingredient we need is an expression for the cutoff value y w (A) \ recall that this marks the y 
boundary below which a given kick will not clear the separatrix U with radius ru(y) = y/l — \fy + 1. 
This quantity does not depend on the slow dynamics and is therefore the same for both our cases. 
For a kick amplitude A, y w (A) = (1 - A 2 ) 2 — 1 if A e [0, 1] (a weak kick). When A > 1, the kick is 
strong and we set y w ysN — 1- 

Using the expressions above we can build our kick map by using ^ and rescaling time by the 
period of a full cycle. To find the period for a given set of parameters, we use h^ 1 and hp 1 to derive 
the silent and active passage times T# and T P . We reiterate that the only dependence on a kick's 
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Figure 6: Analytically (solid lines) and numerically computed (markers) kick maps Fa(0) for three 
values of kick amplitude A. (a) Model parameters: {e = 0.01, w — 1, a — 0.8, b = 0}. (b) Model 
parameters: {e 0.01, w — 1, a — 0.4, b 0.5}. 



amplitude A is implicitly contained in the expression for 9 W = h^ 1 (y w (A)). As mentioned above, 
whenever A > 1, the kicks are strong (6 W = 0) and the resulting maps have a fixed shape. On the 
other hand, as A decreases below 1, the weak kick effect progressively uncovers the left branch of 
the map (9 W ^ 0). 

To verify our derivation of the kick map, we choose 3 prototypical kick amplitudes: one strong 
(A = 1.5) and two weak (A = 0.5 and A = 0.1). We then plot the associated maps for two distinct 
set of model parameters: {e 0.01, w — 1, a = 0.8, 6 = 0} and {e 0.01, w — 1, a = 0.4, b 0.5} 
in order to better visualize the effects of linear (b = 0) and saturating exponential (b > 0) ramps 
for the y dynamics. 

In Fig. |6| we plot numerically and analytically computed maps, rescaled to the unit circle. 
Observe that we get very good agreement between the two and that the main features of weak 
kicks are captured by our phase reduction model. However, numerically computed maps have fine, 
plateau-like segments. As argued in [7J, there are as many of these plateaus as there are spikes in 
a burst. They appear since a kick can induce bursts with spike counts ranging from one to the 
number seen in an unperturbed cycle. Since these numbers are integers and we numerically identify 
phase zero with the last spike of a burst, these plateaus are formed by phases that induce the same 
number of spikes following a kick. This is not captured by our analytical derivation of the kick map, 
which is computed from averaged conditions on the fast variable. However, we will see in Sect. [4] 
that the presence of noise in the system tends to diminish these plateaus, and in Sect. 5 that key 
aspects of synchrony and desynchrony for the original ODEs are predicted by our derived kick map. 

We close this section with some remarks concerning the generality of features found in the kick 
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map. We use the normal form model ([2]) for its analytical tractability but the general mechanism 
responsible for bursting, and hence the associated kick map, have characteristics that span across 
models. For example, for both linear {b = 0) and saturating exponential (b > 0) slow dynamics - 
which imply distinct interactions between the slow and fast subsystems - the maps are qualitatively 
identical. Specifically, at the end of Sect. |2.4[ we noted that region of expansion in the left part of 
the weak kick map is a general feature of elliptic bursters with slow passage effects. The fact that 
the right branch follows the identity is another general attribute. 

The shape of the middle branch (read left branch for strong kick map) is, however, more model 
dependent. Notice that both of our strong kick maps (for A = 1.5) have a left branch that steepens 
as we move leftwards. This is due to the dependence of hp 1 on 7p(y), which grows as y/y. As a 
result, for A sufficiently large, < — 1 for 6 G (0, 9 C ) where 9 C is the root of 



For our parameter set with b = 0, 9 C ~ 0.0968 while for the second set with b > 0, 9 C ~ 0.0619. 
This curvature shrinks as we decrease the parameter a in ^ and the silent phase lengthens. 

In general, this branch depends on the slow dynamics and the ratio of silent to spiking times 
Ts/T P) which impose the following constraint: by construction, a strong kick implies Fa(0) = 1 
and Fa(Ts) — T$. If we were to approximate the y dynamics by constant velocities (as done in [7]), 
this map branch would be linear and hence contractive, whenever T$ > Tp. In general, as long as 
the latter is true and the functions hs, hp (desribing y dynamics) have small total variation, we 
can expect this branch to be mostly contractive, as for both cases explored above. We note that 
we obtain such contraction for the biophysical systems we study in Sect. |5.2[ meant to model GPe 
neurons in the Parkinsonian state; [7j draws an interesting contrast with cases having T# < T P . 
The generality of these features motivate the analysis of dynamics induced by the kick map as we 
show in the next section. 

In light of these remarks and in the interest of clarity, we use the maps computed above for 
the parameter set with b = to carry out our analysis for the rest of the paper. We stress that 
the arguments that will follow hold true for other parameter sets of system ^ and any elliptical 
bursting system having the features described above. 

3 Dynamics of the kick map 

Now that we have an understanding of an elliptic burster's response to input kicks of various 
strengths A, we turn to the other input parameter of relevance - the period between these kicks r 
- and study the iterated dynamics of the map for various combinations of A and r. 

3.1 Iterative framework 

We now use the kick map to build an iterative dynamical system capturing the evolution of cells 
subject to periodic stimulation. Let 



d 
d9 



hp 1 o h s (0 c ) + 1 = 0. 



(22) 



F Ar (9) = F A (9)+r (mod 1) 



(23) 
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which returns the phase of a kicked cell right before the next kick, r time units later. The parameter 
r translates the map vertically as illustrated in Fig. [7| For a chosen pair (A, r) and some initial 
phase O ? an orbit is defined by 9 n+ i = F^ r (# n _i) = FjJ (0q)- 

We note that, in contrast with the map of [7J, our kick map does not become rescaled as r is 
varied, but is rather translated around the circle (as in related studies 0EE]). Another difference 
is the presence of sustained expansion even though T$ > T P (see above), due to accelerated slow 
passage effects induced by weak kicks; we will show that this expansion leads to positive Lyapunov 
exponents for certain values of r. This phenomenon has been exploited in simpler dynamical systems 
with slow passage through a supercritical Hopf bifurcation in the context of chaos control [54] . 




Figure 7: Effect of kick period r( mod 1) on the kick map. Original kick map Fa(0) in black and 
r induced kick map Fa, t (@) m re d. r simply translates the map vertically. 



Equation (23) assumes that a kick acts on cells located on unperturbed trajectories. As also 
noted above, the separation of timescales for elliptic bursters implies very fast attraction back to 
steady state trajectories following a kick, so that this assumption is generally valid. Moreover, the 
iterated dynamics for small values of r are relevant in any case, as they can be accessed by longer, 
equivalent kick periods modulo one. 



We next use Eqn. (23) to study the response of a population of identical elliptic bursters with 
different initial conditions to a common, pulsatile signal. For example, globally stable fixed points 
or periodic orbits represent phase locking regimes towards which the long term behavior of any cell 
will converge. On the other hand, maps that yield sensitivity to initial conditions and complex 
orbits are representative of desynchronizing inputs, when delivered to a population. To further 
explore population behavior, we need to define a metric by which we quantify synchrony of phase 
points on S 1 . We call this our "synchrony measure" and describe it next. 

3.2 Assessing synchrony 

There are many ways one can quantify how closely N points are distributed on S 1 . Two natural 
choices are the binned entropy H and order parameter R (also known as vector strength) : 

1 N 



R{O u ...,0 N ) = \±£e i »* e * 



3=1 
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For iif , we divide S 1 into TV equal length subintervals, or bins, and take pj to be the number of 
phases in bin j over TV (using the convention OlogO = 0). H takes its maximal value one when 
there is a phase spread into each bin and its minimum value zero when all are concentrated into a 
single bin. On the other hand, R takes its minimum value zero when phases are evenly distributed, 
and one when they are all equal. 

Each measure has strengths and weaknesses as a metric of synchrony. For example, R can be 
zero if the phases are split into two equal, antipodal groups on S 1 - this is hardly an asynchronous 
state. For a large AT, H can take relatively large values even if cells are distributed in close by bins. 
By taking 

w(e 1 , e N ) = ^[R{6 1 , ...,e N ) + (i - h(o u e N ))]. 

to be the average of R and 1 — if, we can be assured that low measures of W correspond to 
cases where cells are well distributed across bins and that these bins are broadly spread around S 1 . 
Specifically, we assess synchronization properties of a given map by taking TV cells {9 n }i< n <N with 
some initial distribution on S 1 , pushing these states forward through m iterates, and computing W 
as the average over the last k out of these m iterates: 

w = lY. w{F i (e l ),... 1 F\e N )). (24) 

i=m— k 

Throughout the paper, we take m > 100 and k = 20, having found empirically that values change 
little with larger values of either. 



3.3 High period orbits and positive Lyapunov exponents 

We next investigate how iterations of our kick maps act on a population of cells for the three 
prototypical cases of strong (A = 1.5) and weak kick maps (A =0.5, 0.1), plotted in panel (a) 
of Fig. [6j We note that, while smooth maps on the circle are well characterized [621 |44j, the 
discontinuities in our map introduce a number of distinct phenomena - such as border collision 
bifurcations with period adding at all orders, and "sharp" transitions to chaos. There is an ongoing 
effort to build a theory to better understand such systems [6j El 1341 14T] . 

We first plot orbit diagrams with respect to the parameter r for each of the three maps at hand 
(Fig. [8]). Specifically, we select 100 cells uniformly distributed on S 1 ] for various r G [0,1], we 
compute the positions of these cells after 150 iterates of Fa jT (0) and "vertically" plot the result. 
Directly below these orbit diagrams, we plot the synchrony measure W of these end states, together 
with numerically computed Lyapunov exponents A for each r, averaged over all trajectories. Finally 
we compute approximations of invariant measures for each r, using a variation of Ulam's method 
developed in [22] (this produces a discretized approximation of fixed densities for a map's Perron- 
Frobenius operator). The results are plotted in the bottom panels of Fig. [8j 

We begin by describing results for the strong kick map (A = 1.5). Recall that the leftmost part 



of the strong kick map ,when 9 G (0, C ) (9 C ^0.0968, derived in Eq. (22)) has a derivative greater 
than one in absolute value. The derivative is less then or equal to one in absolute value everywhere 
else. It is easy to see that the map has a fixed point for any r. If the map intersects the identity at 
9 G (0, C ), the fixed point is unstable and we find that complex dynamics emerge. We can witness 
this by looking at the orbit diagram of this map for small r. Nevertheless, we see that cells tend to 
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Figure 8: Orbit diagrams and related measures for kick maps (Eqn. ([9])). Top to bottom: orbit 
diagrams of 100 cells after 150 iterations; synchrony measure W (red) and averaged Lyapunov 
exponent A (blue); invariant measure approximates. Left to right: strong kick with A = 1.5; weak 
kick with A = 0.5; weak kick with A = 0.1. Marked values for tq below which Lemma [T] applies. 



cluster in small regions and hence are relatively synchronized. For any other r ^ 0, we have stable 
fixed points, implying phase locking of cells to the input kicks. Note that for all maps, r = implies 
a continuum of neutrally stable fixed points, as the right branches of the maps align perfectly with 
the identity. 

For the weak kick maps (A — 0.5, 0.1), some values of r yield stable, discrete attractors as 
well: fixed points and periodic orbits. More interesting are values of r which produce thick, chaotic 
attractor-like objects. These are associated with what appears to be locally absolutely continuous 
invariant measures and positive Lyapunov exponents. It is in these regimes that our synchrony 
measures show a greater spread of cells. This is indicative of chaos, the presence of which is 
consistent with expansive regions of weak kick maps. Interestingly, these regions appear wether a 
map has a gap or not, as we show below. 

Rigorously assessing the presence of chaos is, however, not a simple task. This was accomplished 
for related piecewise smooth maps [9J [HI [45J [12] where various definitions of chaos were used, de- 
pending on context. For example, Keener [45] showed that piecewise, surjective and non-decreasing 
maps with overlap have rotation numbers spanning a non-empty interval, an indication of chaos for 
circle maps. Unfortunately, our weak kick maps sometimes fail to have surjectivity (e.g. A=0.5) 
and always fail to be non-decreasing. One can also try to define trapping regions and a family of 
intervals for which interval images cover at least one other interval and the image of at least one 
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interval covers at least two others. This constructs an shift on a space of sequences, which can 
characterize a chaotic system. However, building such an interval family for the weak kick map 
proves to be quite complex as intervals get flipped by decreasing parts of the map and severed 
by discontinuities. Thus, we do not aim at a complete characterization of the complex dynamics 
produced by our weak kick maps. However, we next show that, for some r values, there is a positive 
lower bound on (sustained) Lyapunov exponents of any trajectory. 

In Fig. [8j we see that some values of r (0 < r < 0.2) induce dynamics that appear chaotic for 
both of our weak kick maps. Thus motivated, we state the following lemma: 

Lemma 1 Consider a piecewise defined map F{9) on the circle which is smooth on three non- 
intersecting intervals I\, I 2 and Is with I\ [J I 2 [j Is = S 1 . Suppose l^flij > a > 1, \^\i 2 \ > b > 
- such that ln(a) > |ln(6)| - and additionally that |/ 3 = 1. Then if F(I 2 ) C 7 3; F(Is)f]I 2 = 
and F(I 3 ) ^ Is, the Lyapunov exponent associated with the orbit of almost any initial condition 
9o G S 1 J if well defined, will be strictly greater than zero (X(9 ) > 0). 

Proof: Given 9 and its forward orbit {0 n }n=o,i...5 the local Lyapunov exponent can be written 
as A(0o) — li m iv^oo J2 n =o l n IJj^X^OI- The derivative is defined everywhere in S 1 except the 
4 border points of the intervals /i,^-^ which is obviously a measure zero set. The condition 
F(I 2 ) C Is imply that any point in I 2 is sent to 7 3 . Since F is smooth on 7 3 and ^f|/ 3 = 1, 
the condition F(I 3 ) ^ Is implies that any point in 7 3 must eventually exit it. Let C = max^jn = 
min m {m | F m {9) £ 7 3 } | 9 G Is}- Then any element of 7 3 stays in 7 3 at most C iterates. Furthermore, 
F(Is)f]I 2 = implies that elements of 7 3 are eventually sent to I\. When an orbit point visits 
7i, it contributes at least ln(a) > to the sum in A(# ), at least ln(6) (possibly< 0) if it visits I 2 
and when it passes by 7 3 . By tracking which intervals an orbit visits, any admissible subsequence 
featuring I 2 must contain I 2 —> 7 3 . The sequence that contributes the least to X(9 ) is therefore 
h h ••• h 1 1 where 7 3 is repeated C times. It follows that for almost every Oj 
A(^o) > 2Tc(ln(a) + ln(6) + + ... + 0) > 0. □ 

We note that weaker conditions could be stated under which a positive Lyapunov exponent 
results, but the above are sufficient for the map at hand, as we now show. In particular, we show 
that Lemma [T] can be applied to the weak kick maps for certain values of r. Clearly, the intervals 
(0, 9 W ), {9 W) Ts) and (T#, 1) will play the roles of I\J 2 and 7 3 . For our prototypical weak kick maps, 
we have I 1(0,6^)1 > 2.7, \ dF ^ T \ (o Wi T s )\ > 0.65 and dF ^ r |(r 5 ,i) = 1 which fulfills the derivative 
criteria. It remains to show that the intersection requirements for interval images are met; we take 
a graphical approach which leads to conditions on r. In Fig. [9j panels (b) and (c) show the map for 
A = 0.5 along with the marked intervals for two distinct values of r and cobweb diagrams of sample 
trajectories. One can verify that the Lemma's assumptions are respected as long as the leftmost 
tip of the middle branch stays smaller than one and the rightmost tip is greater than the identity 
(panel (c)). This happens when < r < tq where 

r c = l- lim F Afl (6). (25) 

Although we do not explicitly graph it, the same argument holds for A G (0, 1). We analytically 
compute these upper bounds for r and get tq — 0.205 for A = .05 and tq — 0.26 for A = 0.1. 
As expected, these are close to 0.2, the rough higher bound we predicted earlier from Fig. |5J In 
addition to positive Lyapunov exponents, r <tq imposes a cyclic structure where trajectories visit 
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a large portion of all three intervals in finite time, which is necessary for the population to become 
widely distributed around S 1 . 

For r > tc, various dynamics can be observed. Figure [8] shows that positive Lyapunov expo- 
nents can still be found sporadically but the regions on which the trajectories accumulate can be 
considerably smaller. Periodic orbits of various periods are also present; in particular, high r values 
seem to be associated with stable fixed points. This motivates a separation of the (A, r)-space into 
three regions: I, where Lemma [T] applies; II, where dynamics are complex and transitions between 
what appears to be periodic orbits and smaller chaotic regions can be seen; III, stable fixed points 
(1:1 phase locking). Panel (a) in Fig. [9]shows these regions. We stress that further analysis of region 
II might yield additional structure but the complexity of our map (circular domain, increasing and 
decreasing parts, transitions from gaps to discontinuity points, etc.) renders a complete analysis 
outside of this article's scope. 

From what was described above, we see that region I is given by 

I : {(At)\0 <t< t c {A) = 1 - lim F A ^9)} 

where we have written the expression for tq to highlight its dependence on the kick amplitude A. 
As A increases, W decreases and consequently, tq as well. Region I vanishes altogether when the 
kick transitions from weak to strong - at A = 1, where tq = 0. 

To define the boundary between regions II and III, we need to derive a condition under which 
our maps have a stable fixed point. We begin by inquiring about when the middle branch of the 
map intersects the identity. This happens when Tc + W < r < 1. The stability of the resulting 
fixed point depends on the derivative of the branch at the intersection point. 

Again, as A increases and both 9 W and tq decrease, more of the middle branch of the map 
is exposed. When r = tq + 9 W) the leftmost tip of this branch intersects the identity. If A is 
sufficiently large, the resulting fixed point will be unstable; recall that for the strong kick map 
(A > 1), |^| > 1 when 9 G (0,# c ). It follows that if 9 W < # c , the first fixed points to appear 
as r increases are unstable. We include unstable fixed points in region II and get the following 
definitions 

II :{(A,t)\t c (A) <t< T C (A)+max{e w (A),d c }} 
lll:{(A,T)\T C (A) + max{6 w (A),e c }<T<l}. 

What can be taken from our analysis thus far is that to desynchronize cells with a r-periodic 
input, the best strategy appears to be weak kicks with < r < tq (region I). To illustrate this, 
panels (b) and (c) of Fig. [9] show sample trajectories for 50 cells under the action of the weak kick 
map with A = 0.5 in both synchronizing (region III) and desynchronizing (region I) regimes. To 
mimic a synchronous state, our initial phases are drawn from a uniform distribution on an interval 
of width .01, around 9 = 1/10. The first 50 iterates are taken with respect to the identity map, 
reflecting a preliminary period with no pulsatile inputs. The next iterates are taken with respect 
to the weak kick map with r = 0.5 (region III, panel (b)) or r = 0.1 (region I, panel (c)). 

As expected, for r = 0.5, every cell attracts to a single fixed point and the synchrony measure is 
maximized. However, for r = 0.1, the cells quickly become widely distributed. This was foreseeable 
from the shape of the computed invariant measures at this value of r. Importantly, the synchrony 
measure drops substantially. 

We close this section with an important remark concerning circle maps and phase locking. Here, 
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Figure 9: (a) Three regions of map dynamics in (A,r)-space (see text for details), (b) & (c): two 
sample evolution trajectories for the kick map with A = 0.5. Top to bottom: cobweb diagrams and 
marked intervals (O,0 W ), (6 w ,Ts) and (T#, 1) for which Lemma [l] applies; sample trajectories with 
50 initial condition randomly chosen close to 6 = 0.1 (iterations 50 and above are from kick map); 
synchrony measure W. (b) r = 0.5, cells synchronize; (c) r = 0.1, cells desynchronize. 



our kick map can be seen as a A-perturbation of a r-rotation. Typically, smooth perturbations 
of rotations admit Arnold Tongues: well separated wedge like regions where phase locking with 
rational rotation number occurs (see e.g. [67J, Sect. 21.6). In our case, the loss of smoothness 
introduces sustained region where chaotic dynamics prevail, hereby showing surprising effects of 
discontinuous circle maps. 



4 Effects of noise on kick map and synchrony 

Up to this point, we developed a phase reduction framework to analyze the dynamics of periodically 
forced bursters. We now ask: does the phase reduction remain valid in the presence of noise? If 
so, what qualitative changes in the map occur, and what are the consequences for entrainment of 
bursters? In this section, we develop answers through a blend of numerical results and analytical 
approximations. 

The effects of stochastic perturbation on elliptic bursters have been previously studied in various 
contexts [21 [63], [48]. The unifying theme is the effect of noise on slow passage through a Hopf 
bifurcation. Here, we build on these results to better understand the role of noise on burst responses 
to periodic pulsatile inputs. We will comment further on prior results as we progress. 

To this end, we introduce a stochastic perturbation in the fast variables of the normal form 



21 



model, so that Eqn. ^ becomes 

z = (y + iw)z + 2z\z\ 2 -z\z\ 4 + I(t) + tf(t) 

( I i2 7 \ ( 26 ) 
y = e(a- \z\ -by). 

where 77 > is a small, real, noise strength parameter and £(t) is a time-periodic train of discrete 
small kicks with normally distributed amplitudes. That is, £(t) = — zAt), where the & are 

i.i.d. as AT(0, y/At) and the timestep At controls the temporal resolution of the noisy perturbation. 
We take At = 0.05, the numerical solver's maximal timestep, so that approximates white noise 
(see Sect. [5] for more on numerical methods). Notice that the noise term only acts on the real part 
of the fast variable z, which mimics a cell's voltage variable. 

In [2j , the authors show that such noise terms diminish slow passage effects through supercritical 
Hopf points (i.e., causing cells to jump to the spiking state closer to in [63], a similar effects 
was found for elliptic bursters [63, 48j. As we discuss below, there are cases where a phase reduction 
can still be defined in the presence of such noise, with an interesting and tractable impact on the 
kick map's shape. 



4.1 Effects of noise on the phase reduction of elliptic bursters 

Our previously derived phase reduction relied on an important assumption: periodicity of the burst 
cycle. As discussed in Sect. [2| elliptic bursters do not necessarily have periodic solutions, but rather 
a metastability property which guarantees the constancy of cycle's duration T, up to 0(e). This 
regularity is what enables our phase reduction. 

When noise is added to the fast subsystem, either regular or highly variable burst durations 
can result, depending on noise strength and type j2j [63j [48]. Below, we will show that there is 
a wide range of noise strengths that significantly impact the underlying dynamics, but maintain 



regular burst durations. Specifically, for system (26) (e = 0.01, w = 1, a = 0.8, b = 0) we compute 



the coefficient of variation (CV) of these durations as described in Sect. 2.3, for noise strengths rj 
ranging from 10 -17 to 10 _1 . 



Panel (c) of Fig. 10 shows our findings. Although the mean period (T) decreases with increasing 
noise strength, the CV remains low - below 10 -2 - for the range of noise strengths rj < 10 -3 (recall 
that CV ~ 10 -3 for the noiseless case). In other words, for a wide range of noise strengths, random 
forcing does not introduce substantial variability to the burst period. This is consistent with results 
from [63] where the authors show in a closely related setting that the silent and spiking times T$ 
and Tp are related to the log of the noise amplitude. 

Despite the fact that they preserve regular burst periods, noise strengths rj < 10 -3 have a strong 
impact on the slow passage effect. The integral condition ^ reflects cancellation of attraction to 
S by and repulsion away from it; noise limits the extent of attraction and therefore the duration 



required for repulsion. This is illustrated in panel (a) of Fig. 10, which shows that the jump up 
point from silent to spiking regimes decreases as rj increases within [0, 10 -3 ]. In this regime, the 
averaged y dynamics are relatively unaffected by this stochastic forcing, so that CV remains low. 
Once r] = 10 -2 , there is no slow passage and the solutions jump up to spiking as soon as y crosses 
yn. For r] > 10 -2 , the stochastic kicks have accumulated effects comparable to our forcing kicks and 
we see solutions randomly jumping into spiking before y reaches which explains the increasing 
CV for this range of noise strengths [63j . 

We next pursue phase reduction to an approximate, deterministic circle map for the low CV 
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cases (77 G [0, 10~ 3 ]). 



4.2 Effect of noise on the kick map 



For strong kicks - where responses are not determined by slow passage effects - adding noise does 
not considerably change the shape of the kick map; the chief effect is that the region where |^| > 1 
shrinks (not shown). This similarity was expected since a strong kick instantaneously translates a 
bursting solution to the active phase, where noise has little effect since fast dynamics follow large 
amplitude trajectories. 

For the remainder of this section, we concentrate on the considerable changes noise induces for 



weak kick maps. In panel (b) of Fig. [TO] we plot numerically computed maps for A = 0.5 and r] = 

-9 



10 _ib , 10 _y , 10 . Each marker represents the average response of a given initial phase to 10 distinct 
realizations of stochastic forcing. In what follows, we derive approximations for these maps, plotted 
in solid lines on the same figure. In order to proceed, we discuss key differences among numerical 
maps of stochastic bursters with varying noise strength. 




Figure 10: (a) Slow passage effect shortens with increasing noise, (b) Numerically and analytically 
computed kick maps ( from Eqn. (36)) for noise strengths 77 = 10 -15 , rj = 10 -9 and rj = 10 -3 (all 
with kick amplitude A = 0.5). (c) CV of burst cycle period (computed using 150 cycles) against 77 
using a log-log scale, (d) Average burst cycle period (T) (same sample as in (c)) against 77 using a 
log-linear scale. 



We can see that the maps are qualitatively similar in almost all aspects except the left branch, 
which collapses on the identity as noise increases. Recall that the expansion in this branch is due 
to altered slow passages due to weak perturbations on the stable part of the silent branch S. As 
discussed above, noise shortens slow passages and one might expect the steepness of the expansive 
part to decrease as noise increases. Surprisingly, it is not this steepness that changes but rather 
the onset of the expansive ramp {0b defined below) which varies. To explain this phenomenon, we 
must turn to the concept of buffer points for delayed bifurcations. 

When deriving the analytical kick map ^ , we relied on the integral condition ^ which dictates 
that the further a solution starts (at yi) from the Hopf point the longer the slow passage will 
be. This is true for yi in a range before up to a point beyond which this relation fails. In 
fact, given some y dynamics, one observes that the length of the slow passage (yj — y#) saturates 
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to a constant value for any initial yi far enough from y H . This is called the maximal delay. The 
closest point to the "static" bifurcation point (yn) to generate a maximal delay is called the buffer 
point of a delayed bifurcation [21J. The techniques used in [21 [63], HHJ |3] to analyze deterministic 
and stochastic slow passages assume that a solution always remains closer to the bifurcation than 
the associated buffer point. This is also the case for the bursters we study in the absence of noise, 
where the buffer point is to the left of ysN- Thus, it does not affect the dynamics, enabling us to 
use condition 



For noisy bursters, the numerically derived weak kick maps in Fig. [10| show that the buffer point 
can lie to the right of ysN, therefore exerting an important effect. If a weak kick is delivered when 
yi is to the left of the buffer point, the trajectory retains no memory of the kick and the phase of 
the cell is left unchanged by the kick. If it is delivered when yi is to the right of the buffer point, 
the kick will shorten the slow passage, as for the noiseless case. 

To our knowledge, there are no general results in the literature deriving buffer points for noisy 
delay bifurcations. In what follows, we derive an approximation for them, in the context of weakly 



kicked trajectories, that holds for Eqn. (26). We work in the small At limit, so that the stochastic 
forcing is white noise. 

4.2.1 Buffer points and stochastic slow passage through the Hopf point 

For a given noise strength 77, define ys to be the buffer point for the delayed Hopf bifurcation in 



Eqn. (26). That is, y# is the smallest initial value yi such that a weak kick delivered at yi will 



induce a change in the jump up point yj. Naturally, y# < yH and judging by the shape of our weak 



kick maps in Fig. [TOj we can expect ysN < Vb- We now derive a probabilistic criterion that gives 
an accurate approximation for y B . 

We begin by assuming that At is sufficiently small in Eqn. ( [26] ) so that the noise term rj£(t) can 
be approximated by the white noise term rjdW(t)/dt, where W(t) is a real valued Wiener process. 
We are interested in solutions following the silent branch S when y 6 (ysN,yii)- As done in Sect. 



2.5, we make the singular limit assumption that z travels on S and therefore that \z\ = 0. This 



enables us to decouple the y dynamics and get an expression for y(t) = yi + g(st), precisely as in the 



deterministic case (see Eqns. (16),(|17|)). These assumptions hold in the limit that s and 77 — ^ 0, and 
we will show that they give good approximations for the parameters used here (e = 0.01, 77 < 10 -3 ). 

Since we are interested in solutions near S, which is composed of fixed points of the fast subsys- 
tem, we linearize the fast dynamics about z — 0. We write the resulting equation in real coordinates 
x = (xi, x 2 ) T where z = X\ + ix 2 : 

dx = J{t)xdt + BdW{t) (27) 



where 



J{t) ~ ' w y(t) h B ~ I 



and W(t) is now a two dimensional real valued Wiener process. Notice that J(t) commutes with 
itself (J(t)J(s) = J(s)J(t) V £, s) which enables us to write the noise-free (77 = 0) solution of (27) 

as 

x(t) = eti J ( s)ds Xi 

where x(0) = X{. Using this property, Ito's formula yields explicit expressions for the mean ji{t) 
and covariance matrix E(t) of the time dependent probability distribution for x governed by Eqn. 
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(27): 



fi(t) = efi'^fjLi (28) 

pt 

E(t) = eti J{s)ds Zieti J{s)Tds + / dt'efi J{s)ds BB T eti J{s)Tds (29) 

Jo 

where /i(0) = and E(0) = E^. See [26], Chap. 4 for details of this derivation. In particular, we 
will suppose that we have an initial distribution for x that is Gaussian; then, the distribution of 



x at any time t is fully determined by Eqns. (28) and (29), as a Gaussian distribution with mean 



n(t) and E(t). We introduce the following notation p(-) for this distribution, which makes clear 
the dependence on the initial condition for y(0) = yi as well as the initial distribution of x and the 
elapsed time t: 

x ~ p(x(t)\y h /z»,E») . (30) 

We next study the distribution of our (linearized) fast variable x at the Hopf point yn when no 
kick is delivered. This will give us a reference unkicked, or "natural" distribution p n {x) important 
in computing the buffer point below. The trajectories of interest jump down from spiking when 



y = yi = ysN and take g 1 (yH — Vsn)/£ time units to reach y H . We use Eqn. (30) to write the 
resulting distribution as 

p n {x) = p(x(g~ 1 (y H - y S N)/e)\ysN, (0, 0) T , E SN ) (31) 

where 

v _ ( r P (y SN ) 2 
hsN ~ { r P {y SN f 

and we have used Eqn. to substitute in for time in Eqn. (31). The natural distribution is 
therefore defined to have an initial variance equal to the squared radius of the periodic orbits on P 
when they vanish at ysN, marking the end of the spiking phase. Due to the long timescale of the 
slow dynamics, however, we find that the choice of initial variance has little effect on p n . We note 
that p n {x) is centered at x — (0, 0) T with covariance depending on rj. 

Next, we ask whether a trajectory that has received a weak kick can be expected to undergo 
comparable slow passage through the Hopf point as for the unkicked trajectories described by p n (x). 
Note here that an unkicked trajectory admits a maximal delay going through the Hopf bifurcation, 
and that some kicked trajectories will also have the same slow passage, as we expect ysN < yB- 

If kicked trajectories typically pass through y H at locations x with high probability density 
according to p n (x), we expect that they will have a comparable slow passage times as unkicked 
trajectories. On the other hand, if these trajectories are typically found where p n (x) is low, this 
indicates that they are further away from the branch of equilibria S. Thus, they will tend to escape 
S (i.e., jump up) sooner than the unkicked solutions. We obtain an approximation for the buffer 
point ys by asking when this distinction between kicked and unkicked trajectories at yu occurs. 

To this end, we compute the "kicked" distributions of trajectories Pk,y(%) for which a weak 
kick of amplitude A is applied when y(t) = y. Upon their arrival at yH, we approximate these 
distributions as 

Pk9V (x) = p(x(g-\y H - y)/e)\y, (A, 0) T , E n ) (32) 

where the mean trajectory is at {A ) 0) T following the kick and E n is the covariance matrix fromp n (x), 
under the assumption that the trajectory followed the natural burst cycle before the translation 
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Figure 11: (a) Plot of d(y) with respect to y for rj = 10 -15 , 10 -9 , 10 -3 . Dashed black line marks the 
threshold d# determining the points y# at the intersections with the d(y) curve, (b) Plot of buffer 
phase point 6b with respect to noise strength 77, marking the onset of expansion in the weak kick 
map. Kick strength is A — 0.5, as in Fig. [10, 



induced by the kick. We now assess to what extent p^ y and the natural distribution p n overlap. 
We use the symmetrized Kullback-Leibler Divergence between the two distributions 

d (v) = \ ( D KL[Pn(x)\\Pk,y( X )} + D KL\Pk,y{x)\\p n (x)]) (33) 

where 

Dkl\p\\q]= f p(x)\n^dx. (34) 
Jr2 q{x) 

For y sufficiently far from y H) the distribution p^ y has enough time to converge close to p n (x) 
before y(t) reaches As a consequence, d(y) will be close to zero. We define the A-dependent 
buffer point ys to be the first value of y for which d(y) grows beyond a threshold Since the 
distributions pk, y and p n are sharply peaked Gaussians (with variance of order ?y 2 ), d(y) quickly 
explodes - to several orders of magnitude above one - when the two distributions fail to overlap. 
Therefore, we choose ds = 10 1 as a good indication of separation among kicked vs. unkicked 



trajectories, as illustrated in Fig. 11 (a). 

Finally, we determine the phase point 9b corresponding to ys , via Q: 

e B = £M (35) 

where T is the mean period of the unperturbed trajectory. This point marks the onset of expan- 
sion for the associated kick map. We find excellent agreement of this prediction with numerically 



computed kick maps, as seen from Figs. 11 (b) and 10 (b) 



4.2.2 Deriving the kick map for noisy bursters 

We next derive an approximate expression for the complete kick map in the presence of weak noise. 
Our first step is to decouple the y dynamics from the (noisy) fast variables. This is guided by the 
assumption that, due to weak noise, most trajectories closely follow S and P. We then proceed to 
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derive hs, hg 1 and hp 1 as done in Sect. 2.4, with a single modified value: yj. Indeed, as described 
above, the jump up point yj is closer to yn for noisy bursters and we numerically compute its value 
for each noise strength 77. 

In the presence of a buffer point, a weak kick can now have two outcomes: either it has no effect 
if it is received when < y < y# as it does not alter slow passage, or it shortens slow passage as 
in the deterministic case, when ys < y < y w - To capture the slow passage effects of weak kicks 
(responsible for the expansion in our map) in the presence of noise, we still use integral condition ([8]). 
Although this formula was derived for deterministic bursters, it relies in the linearization of the fast 
dynamics about S : which we assume remains valid even in the noisy case. As a result, numerical 
simulations show that yj(y) holds true for y G (vb^Vh) except in a short interval to the right of 
ys where small errors are observed. Note that these errors diminish with smaller noise for which 
y^ys) is quite close to the numerically computed yj. We proceed to write an expression for our 
new kick map, which now contains an addition piecewise-defined section arising from the presence 
of 6 B : 

e if e e [o, e B ] 

F(0 \- \ ° + h ~p 1 (yj(hs(e)))-h- s 1 (y j (h s (e))) aee[e B ,e w ) 

fc-ioMfi) i£6e[d w ,T s ] (36) 

if0e[T s ,l]. 



We obtain excellent fits as shown in Fig. 10 (b). We end by noting that we get similar fits for 
various kick strengths A as well as distinct parameters sets (i.e. b > 0) for the normal form model 
(not shown). 

4.3 Effect of noise on iterated dynamics 

We now explore the dynamical properties of the kick maps computed in the presence of noise. 
Figure p~2] shows orbit diagrams, synchrony measures and averaged Lyapunov exponents three maps 
(computed as for Fig. [8]). For r\ = 10 -5 , 10 -9 , the maps retain some expansion and we see behavior 
that appears chaotic for small positive values of r. This range of r values shrinks as the as the 
expansive region of the map gives way to neutrality with increasing noise; at the same time, the 
"support" of the orbit diagrams appears to decrease. Arguments similar Lemma[T]can be formulated 
for these cases to show the existence of positive Lyapunov exponents (see also bottom panels of 
Fig.§. 

For the map with rj = 10 -3 , there is no expansive region: although there is still a slow passage 
effect the kick's amplitude A = 0.5 induces a cutoff 9 W small enough such that the system does not 
retain memory of any kick (9 W < 9b)- Thus, we cannot expect positive Lyapunov exponents (see 



bottom panel in Fig. 12). However, the orbit diagram shows a broad spread of points for small 
positive r. These are stable, high period orbits, originating from border crossing bifurcations as r 
increases. This has been established by Belair in the context of periodically forced integrate and fire 
oscillators [5], where a very similar map is studied: he shows that as the map is shifted vertically, 
stable periodic orbits of a wide range of periods can be found, following a Farey tree sequence. 

In sum, small positive values of r result in either positive Lyapunov exponents or high period 
orbits for the weak kick maps, depending on the noise level rj assumed in deriving the (deterministic) 
map. In the first case, expansion directly desynchronizes cells; in the second, we will see that the 
high period of orbits, coupled with additional variability due to the underlying noise, can have a 
similar effect. 

We now introduce stochastic terms into our discrete kick map dynamics to account for the 
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Figure 12: Orbit diagrams and related measures for weak kick maps with A = 0.5, for various 



noise strengths r] (see Fig. 10). Top to bottom: orbit diagrams of 100 cells after 150 iterations; 
synchrony measure W (red) and averaged Lyapunov exponent A (blue). Left to right: low noise 
strength r\ = 10 -15 ; medium noise strength r) = 10 -9 ; high noise strength rj = 10 -3 . 



variability in burst periods discussed above (i.e., CV ^ 0). If, for a given cycle, a cell has a 
shorter/longer period than the one used to compute it's kick map, its phase following a kick will 
be slightly shifted from the phase given by the iteration of the map. To capture this, we introduce 
jitters', additive stochastic terms acting on r, independent for every cell. The goal is not to capture 
the exact phase response of cells, but rather to give a qualitative account for the impact of period 
variability on statistical metrics such as our synchrony measure. We proceed as follows: 

Since the construction of our kick maps rescales the period of any burster to be 1, the CV can 
be interpreted as the standard deviation of a burst cycle's period. We define a jitter ( to be random 
variable drawn from a normal distribution with zero mean and standard deviation equal to the CV 
of the case we are considering. Suppose we want to model the phase evolution of M cells subject 
to a common, periodic kick train of period r(mod 1), of amplitude A, subject to stochastic forcing 
of strength rj. The phase of the m th cell right before the n + 1 st kick can be written as follows 

C+i = F AMr+&(^n) ( 37 ) 

where £™ ~ ud AT(0, CV V ). In other words, at every iteration, we draw a different jitter ( for every 
cell. Note that we modified our notation to emphasize the map's dependence on noise strength r] 



(Fig. 10). 



Using (37), we again iterate 100 cells 150 times with added jitters and plot the orbit diagrams 



and synchrony measures in Fig. 13, for the same three levels of rj as in the preceding figures. 
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Jitters, as expected, "smear" orbit diagrams, with a greater effect for larger rj. In particular, note 



V=0, CV=0.0054 r] = 10 9 ; CV=0.008 rj = 10~ 3 , CV =0.015 




Figure 13: Orbit diagrams and synchrony measure for weak kick maps with A = 0.5, for various 
noise strengths rj and added jitters ((CV). Top to bottom: orbit diagrams of 100 cells after 150 
iterations; synchrony measure W. Left to right: low noise strength rj = 10 -15 ; medium noise 
strength rj = 10 -9 ; high noise strength r] = 10 -3 . 

the smoothing of periodic points for small positive r in the high noise case (77 = 10 -3 ). Interestingly, 
the twin effects of noise in reducing expansion but increasing cell-to-cell jitter result in comparable 
levels of the synchrony measure W across the three cases. 

We reiterate our main conclusion: although the underlying mechanisms differ across a wide 
range of noise strengths 77, pulsatile inputs in the "weak" kick regime - with an input frequency 
slightly slower than (a multiple of) cells' intrinsic frequencies - will result in desynchrony among a 
population of recipient cells. 

5 Validity of phase reduction, O.D.E. simulations, and a 
neurobiological model 

In this section, we explore the validity of our phase reductions in both deterministic and noisy cases 
- and our analysis of the discrete kick map that follows - by numerically integrating the O.D.E.s 
themselves. Rather than demonstrating a complete correspondence between the kick map and 
solutions of the differential equations, we seek to verify that an informed choice of kick amplitude 
and period, based on the kick maps, does indeed yield the predicted (de)synchrony behavior among 
solutions to the O.D.E.s. Specifically, we show that small, positive values of r (i.e., r G (0, tc) 
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in the deterministic case) with amplitudes in the "weak" regime lead to the greatest desynchrony; 
conversely, large values of r synchronize cells. 

We first consider the normal form system ([2]), and then turn to a biologically detailed neu- 
ronal model for which our main findings persist. All numerical computations were carried out in 
MATLAB. We use the stiff solver ode 15s with both absolute and relative tolerances set to 10 -6 to 
integrate all differential equations; input kicks and additive noise are treated as non- autonomous 
terms by the solver. 



5.1 Normal form model 

Here, we numerically integrate a population of A = 30 uncoupled cells governed by system ^ or 
its stochastic counterpart Eqn. (26), taking the "large" noise value rj = 10 -3 studied above (with 
{e = 0.01, w — 1, a — 0.8, b = 0}). We concentrate on one weak kick amplitude, A = 0.5. In 
each case, we implement periodic kicks that correspond to r = 0.1 and r = 0.8 for the kick map, 
to illustrate desynchronizing and synchronizing behavior respectively. More precisely, we use a kick 
period equal to (1 + r) x T where T is the natural period of the O.D.E.'s burst cycle. This allows 
trajectories at least one natural period to relax toward the unperturbed cycle in between kicks 
(similar dynamics occur for periods T x (n + r), n G N). For r\ = 0, T ~ 465 while for r] = 10 -3 , 
T-337 (Fig.[l0|(d)). 

Figure [14] displays the results via raster plots: for each cell, a dot is placed at the moment 
that spiking terminates (corresponding to phase 9 = 0). We also plot a synchrony measure for the 
simulated population. This is done by assigning phases to each cell relative to their most recent 
spike termination event, as a fraction of elapsed time partitioned in bins of length T. To better 
illustrate the desynchronizing effect of weak kicks with r = 0.1, initial conditions are chosen at 
random with phases at most 2% apart (i.e., an initially synchronized population); to illustrate the 
synchronizing effect of kicks with r = 0.8, initial phases are allowed to be more sparse. In all cases, 
we let the cells evolve without inputs for a few burst cycles, and then begin to apply the pulsatile 
inputs. 

The results agree well with predictions from the kick map: a weak kick of A = 0.5 administered 
at T x 1.1 successfully spreads cells apart while the same kick with period T x 1.8 synchronizes the 
population. Note that we chose these values of r only as informed guesses; other nearby values can 
achieve similar results in both synchronizing and desynchronizing the population. Moreover, results 
for various other kick amplitudes also agree well with the behavior predicted from the associated 
kick maps (not shown). 



5.2 GPe bursting neuron 

We now investigate whether the mechanisms described above will persist for a more biologically 
detailed model. Specifically, we study a 5-dimensional, Hodgkin-Huxley-type model of a neuron 
from the GPe basal ganglia nucleus |66j [7j. This model produces elliptic bursting where the onset 
of spiking is due to a subcritical Hopf bifurcation and a burst termination is due to a saddle node on 
an invariant circle, in agreement with the normal form system ([2]). In detail, the fast variables are 
the voltage V potassium current gating variable n, sodium current gating variable h ) and calcium 
T-current gating variable r. The slow variable is calcium concentration Ca. The equations are as 
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Figure 14: A population of 30 numerically integrated solutions of Eqn. ^ or (26), all receiving 
a common, periodic weak kick input (A = 0.5). Left column: kick period T x 1.1 (equivalent 
to r = 0.1) results in population desynchrony. Right column: kick period T x 1.8 (equivalent 
to r = 0.8) synchronizes the population. Top row: no noise (77 = 0). Bottom row: high noise 
(77 = 10 -3 ). Black dots give the raster plot (see text); red curves plot the synchrony measure W vs. 
time (see text). 



follows: 
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~dt 
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(38) 



-e[I C a + It + k Ca Ca) 



where the / terms represent membrane currents and are functions of the gating variables and the 
voltage; all definitions and parameter values are as in [66J. Additionally, the terms I(t) and rj£(t) 
represent the pulsatile inputs and the noise term, entering as currents. 

Figure [l5|(a) shows that the shape of action potentials and the timescales differ from those of 
the normal form model ([2]); nevertheless, the dynamics have a very similar structure. In particular, 
panel (c) shows a (projected) 3-dimensional plot of a bursting trajectory together with the skewed 
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Figure 15: Plots from model (38). (a) Voltage trace, (b) Calcium trace, (c) Bursting solution in n, 



V, Ca plane in red and separatrix U in blue, (d) Two numerically computed kick maps for model 



(38) with high noise strength rj = 10 3 ; strong kick (A = 30) in black and weak kick (A = 3) in red. 



separatrix £7, computed using the MATCONT package [20]. Figure 15(d) shows two numerically 
computed kick maps for both strong and weak kick amplitudes (A = 30 and A — 3, respectively). 
These maps were computed in the presence of noise with r\ = 10 -3 , which is the largest noise 
strength that keeps CV at O(10 -2 ) for the simulations; maps represent the average phase response 
taken over ten runs with different noise realizations. 

Overall, the structure of these maps is more complex than for the normal form model. In 
particular, "small" plateaus and associated discontinuities are promienent. As for the normal form 
model, there are as many plateaus as there are spikes in an unperturbed burst, and kicked solutions 
that elicit a certain number of spikes in the subsequent burst accumulate in each plateau. However, 
for this model, the slow variable (calcium concentration) varies more during a spike and creates 
bigger gaps between plateaus. As a result, even for a strong kick, certain values of r yield localized 
stable periodic orbits, as opposed to only fixed points. These appear via border collision bifurcations 
due to discontinuities between plateaus (not shown). However, the small amplitude of these periodic 
orbits keeps the cells attracted to them quite synchronized. 

Additionally, the shape of the left part of the weak kick map is also quite distorted compared with 
maps derived from the normal form model. In particular, notice that there are large discontinuities 
close to zero. This is due to the skewed cone shape of the separatrix U : since the neuron model 
does not have the same symmetry as our normal form system, when the solution drops down from 
spiking, it spirals towards the resting branch and some lobes of this spiral come very close to the 
separatrix. When the solution is kicked, even weakly, on the upper part of a lobe, it passes the 
separatrix and jumps to the spiking state; the same weak kick will not have this effect if it is 
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Figure 16: A population of 30 numerically integrated solutions of Eqn. (38) with noise strength 
7] = 10 -3 , receiving a common periodic weak kick input (A = 3). Left column: T x 1.1 (equivalent 
to r = 0.1) results in population desynchrony. Right column: kick period T x 1.8 (equivalent to 
r = 0.8) synchronizes the population. Other plotting details also as for Fig. 14, 



delivered only moments later. The resulting large gaps in the weak kick map add to the complexity 
of the dynamics for low r. 

Apart from these differences, the prominent features observed in the kick map for the normal 
form model remain: the neutral/expanding left branch for the weak kick map, and the contracting 
middle branch and neutral right branch. We repeat the numerical experiment described in Sect. 



5.1, this time only for the noisy case (rj = 10 3 ), and plot the results in Fig. 16, We see the expected 



synchronization and desynchronization from weak periodic kicks (A = 3) with periods equivalent 
to r = 0.8 and r = 0.1 respectively. While cells do not appear to become as desynchronized for 
the r = 0.1 case as in the normal form model, it is reasonable to believe that a more detailed 
analysis of the kick map for the neural model could identify {A ) r) combinations that would further 
desynchronize cells. 



6 Composition of multiple periodic inputs, and an applica- 
tion to DBS 

Above, we showed how weak, periodic inputs can lead to desynchronization for populations of 
uncoupled bursting cells. But how well can such inputs compete with other, synchronizing effects? 
The answer is important in varied applications. A prominent one is Deep Brain Stimulation (DBS) 
therapy for Parkinson's disease. Here, pathologically high levels of synchrony occur among bursting 
cells in the basal ganglia. Synchrony in some basal ganglia areas is in large part driven by common, 
periodic inputs from other areas (see [60j and references therein). A DBS electrode delivers pulsatile 
electrical signals that are designed to mitigate the effects of this synchrony. Thus, two common 
periodic inputs are received by bursting neurons, possibly with competing effect. 

Using the normal form model ([2]), we undertake a brief demonstration of how our results could 
be applied to this setting, for the GPe basal ganglia nucleus that contains neurons believed to be 
elliptically bursting. We do not attempt detailed, biologically complete modeling or aim for direct 
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clinical relevance, and as such note several limitations. First, the source of intrinsic entrainment 
here is a purely common, periodic drive; lateral connections between bursting cells, believed to 
be sparse and weak in Parkinsonian regime [601 [56] . are neglected. Second, GPe is not the most 
common target for DBS in practice, though it has been the focus of several emerging studies |42ll56] . 
Nevertheless, a better theoretical understanding of the interactions between intrinsic and applied 
inputs to the GPe could, in the long term, contribute to the design of signals that desynchronize 
bursting neurons by targeting key instabilities (cf. J33J [321 [19l [25j ) , taking inspiration from similar 
findings for oscillatory neurons [19j ESI SHI EH [49] . 

We suppose that a population of bursters receives a first sequence of synchronizing periodic 
impulses with period T\ and "strong" amplitude A\. The action of these inputs on burst phases is 
given by the kick map Fa 1jT1 (0). As throughout our paper, this returns the phase of a cell following 
a kick, T\ time units later. Aiming to counteract the synchrony due to the first kick sequence, we 
introduce a second series of kicks of strength A 2 . We assume that these have the same period, 
but are delayed by an amount t 2 . That is, the cell receives a ^-kick t 2 time units following each 
yli-kick. We wish to write the kick map that captures the effect of such doublets of kicks. 

In this context, the shift-time following a ^i-kick must be taken to represent the phase of cells 
right before the ^4 2 -kick and the first application of the map must be Fa 1iV2 (0). Similarly, we must 
shift the A 2 map by t\ — t 2 to retrieve the phase of a cell before the next yli-kick. Note that neither 
t 2 nor t\ — r 2 should be too small for this map to be valid, specifically in the presence of weak kicks 
when the cell must have time to enter its spiking phase before the following kick, for the map we 
derive to remain valid. When this restriction is satisfied, the doublet map is given by 




Figure 17: (a) Orbit diagram and synchrony measure for Fa 1 a 2 ,t 1 t 2 while varying r 2 . (b) Top 
: cobweb diagram of Fa 1 a 2 ,t 1 t 2 with r 2 indicated by the red arrow in (a), (b) Bottom : O.D.E. 
simulation of 20 cells initially synchronized by a strong input and then desynchronized by competing 
weak kicks (starting at red arrows). 

An example is shown in Fig. [TTj Suppose we start with an entraining input of strong kicks with 
A\ — 1.5 and T\ = 0.4. We seek to oppose this synchronizing effect with weak kicks of amplitude 



34 



A 2 = 0.5. We use the two first maps of panel (a) of Fig. [6] to build the resulting doublet map 
Fa 1 a 2 ,t 1 t 2 - I n panel (a) of Fig. 17, we compute the orbit diagram of this map (as done in Sect. [3]) 
while treating r 2 as our variable parameter. Using this diagram, we select r 2 = .375 (marked by a 
red arrow), associated with a low synchrony measure. We plot a cobweb diagram of a sample orbit 
in the top of panel (b) of Fig. 17, which clearly demonstrates the destabilizing effects of expansive 
regions in the doublet map. 

We then verify the properties desynchronization predicted by the doublet map by numerically 
solving the underlying O.D.E.s ([2]), for twenty model cells. We begin with initial conditions such 
that the phases are desynchronized and apply strong kicks (A\ = 1.5) at 1.4 x T, where T is the 
natural period of the bursters (i.e., corresponding to T\ = 0.4). As predicted, the cells synchronize in 



response; see the binned synchrony measure rising up to one in the bottom of panel (b) of Fig. 17 



or the raster plots above. After synchrony has developed (red arrows in panel (b)), we "switch 
on" the sequence of weak kicks, leaving in place the original strong kick sequence. Weak kicks are 
applied 0.375 x T time units following each strong kick (i.e., r 2 = 0.375). The desynchronizing 
impact predicted by the doublet map is clear in both the scatter in raster plots and in the drop in 
the synchrony metric W that develops after the weak input begins to be applied. 



7 Conclusion 

We study the behavior of a population of identical elliptic bursters receiving a periodic sequence 
of pulsatile inputs, or kicks. Our aim is to understand which input sequences will result in desyn- 
chronized vs synchronized bursts across the population. Following and extending the approach in 
[7j, we first conduct a phase reduction of the burst dynamics to a circle map, using a slow/fast de- 
composition. This "kick map" depends on two parameters - the kicks' amplitude A and (relative) 
period r - and maps phases from their states just before one input pulse to their states just before 
the next pulse arrives. We next study the effect of varying A and r using a normal form model for 
elliptic bursting (Eqn. ([2])). 

We find that for strong kicks - i.e., with A sufficiently large so that the cell will always be spiking 
following an input - almost any choice of kick period r resulted in 1 : 1 phase locking, and hence 
synchrony across the population. For weaker kicks, we find a rich dynamical structure. In particular, 
the interaction of a weak perturbation with the slow passage effect through a subcritical Hopf point 
induces an expansive region in the kick map. By varying the kick period, we witness the appearance 
of stable fixed points, periodic orbits and regimes with positive Lyapunov exponent. As expected, 
this leads to desynchronization of the population. Overall, we divide the (A, r) parameter space into 
the three regions shown in Fig.|9](a), corresponding to unstable, desynchronizing dynamics, 1:1 phase 
locking, and intermediate, complex behavior, The former, desynchronizing regime is associated with 
relatively weak kicks of periods slightly slower than the natural burst period (0 < t < tc, see 



Eqn. ([25])). 

We also study the effect of stochastic perturbation via noise terms. We find that the phase 
reduction retains its validity but the kick map changes shape, presenting less expansion as the 
noise increased. Importantly, population desynchrony still results from weak kicks with comparable 
values of r in this case, but through a different mechanism than the instabilities that occur for the 
noise free case. Here, desynchrony follows from a combination of high-period orbits and the noise 
itself. Overall, this phenomenon is related to the discontinuous nature of the circle map at hand; 
1:1 phase locking rather than the complex dynamics observed would be expected for small r for 
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many smooth maps [28j [27] . 

We then test the predictions of the reduced circle maps via numerical simulation of the original 
O.D.E. system, finding qualitative agreement. Additionally, we simulate a more biologically realistic 
model of a GPe neuron, and continue to find agreement with the general predictions of our maps. 
Finally, we show that it is possible to use the kick map framework to study the effect of multiple 
sequences of inputs to a cell population. We build an example showing that carefully timed weak 
kicks can compete with an entraining strong input to successfully desynchronize a population of 
bursting cells. 

As a closing remark, we note that the kick map studied here can also capture the effect of 
pulsatile input signals that are neither periodic, nor have a fixed kick amplitude. For any given 
sequence {A n ,r m }, where A n is the amplitude of the n th kick and r n is the delay between kicks 
n and n + 1, the relevant system is the composition of the maps F^ n/7 - n (#). This gives rise to an 
iterated function system (IFS) acting on S 1 . There is a growing body of literature dealing with these 
objects and their application to this problem could eventually help us to understand the behavior 
of bursting cells under arbitrary - and possibly stochastic - stimulation patterns. 
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